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Abstract 

With  the  aim  of  dynamic  simulation  and  control,  a  cell-level  nonlinear  state-space  dynamic  model  of  solid  oxide  fuel  cell  (SOFC)  based 
on  physical  principles  is  built.  It  demonstrates  that  reactant  diffusion  processes  from  gas  flow  bulks  to  triple  phase  boundaries  (tpb)  play  an 
important  role  in  the  dynamic  behaviors  of  SOFC.  The  simulation  shows  that  under  certain  condition,  the  effect  of  the  double  layer  capacitance 
may  be  neglected.  The  phenomenon  of  slow  rise  of  voltage  in  current  interrupt  experiment  may  be  explained  by  the  dynamic  behaviors  of  partial 
pressures  in  the  vicinity  immediately  above  tpbs.  A  new  equivalent  circuit  is  proposed.  The  results  are  demonstrated  through  simulations. 

©  2005  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Solid  oxide  fuel  cell  (SOFC)  is  a  type  of  fuel  cells 
identified  as  the  likely  fuel  cell  technologies  that  will 
capture  the  most  significant  market  in  the  future.  SOFC 
uses  a  special  solid  oxide  cermet  (mostly  yttria- stabilized 
zirconia,  YSZ)  as  the  electrolyte  and  usually  works  at  a 
high  temperature,  in  the  range  of  800-1000  °C  to  reach 
the  electrolytes  ionic  conductivity  requirement.  A  large 
amount  of  literature  has  been  published  about  SOFCs  as 
well  as  other  fuel  cells.  Most  of  the  work  has  focused  on 
investigating  static  electrochemical  characteristics,  such  as 
reaction  mechanisms,  state-of-the-art  cell  components,  new 
materials,  etc.  For  the  purpose  of  dynamic  simulation  and 
control,  the  dynamic  characteristics  of  fuel  cell  must  be 
understood.  In  this  paper,  a  dynamic  model  of  SOFC  at  cell- 
level  is  proposed  to  investigate  the  dynamic  properties  of  fuel 
cells. 

In  order  to  investigate  dynamic  characteristics  of  SOFC, 
some  dynamic  modeling  and  simulations  have  been  per¬ 
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formed  on  molecular  scale  [1,2].  However,  such  micro¬ 
dynamic  models  are  not  suitable  for  macro-applications. 

Early  macro-dynamic  modeling  of  SOFC  was  performed 
by  Achenbach  [3,4].  He  examined  the  transient  cell  voltage 
performance  of  a  cross-flow  planar  SOFC  cell  due  to  the 
temperature  changes  and  the  perturbations  in  current  density. 

The  dynamic  model  given  by  Padulles  et  al.  [5]  included 
the  species  dynamics  on  stack-level  the  first  time.  Zhu  and 
Tomsovic  [6]  adopted  the  model  of  Padulles  et  al.  [5]  for  ana¬ 
lyzing  the  load-following  performance  of  microturbines  and 
fuel  cells.  Sedghisigarchi  and  Feliachi  [7]  combined  Achen- 
bachs  heat  transfer  dynamics  and  Padulles’  species  dynam¬ 
ics,  and  performed  dynamic  simulations.  However,  in  these 
models,  only  lumped  dynamic  behavior  of  species  along  the 
fuel/air  channel  is  considered.  The  processes  of  species  trans¬ 
port  from  flow  bulk  to  triple  phase  boundary  (tpb)  have  not 
been  considered. 

The  model  proposed  in  this  paper  is  a  cell-level  species 
dynamic  model.  It  lies  between  micro-  and  macro-scale.  Be¬ 
cause  behavior  of  stack  is  determined  by  that  of  cells,  the 
cell-level  model  is  a  building  block  for  a  stack-level  model. 
Dynamic  behaviors  of  voltage,  current,  gas  consumption  rates 
controlled  by  different  load  and  partial  pressures  are  demon¬ 
strated  through  simulations. 
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Nomenclature 

au2 

activity  of  hydrogen 

au2o 

activity  of  water  vapor 

ao2 

activity  of  oxygen 

A 

area  (m2) 

C 

concentration  (mol  m-3) 

Cct 

charge  transfer  capacity  (F) 

D 

diffusion  coefficient  (m2  s-1) 

E 

voltage  (V) 

^act,a 

anode  reaction  activation  energy  (J  mol-1) 

^act,c 

cathode  reaction  activation  energy  (J  mol-1) 

F 

Faradays  constant  (=96487  C  mol-1) 

Agf 

Gibbs  free  energy  released  (J  mol-1) 

i 

cell  current  (A) 

h 

current  limit  (A) 

io 

exchange  current  (A) 

js 

consumption  rate  flowing  into  the  outer  surface 
of  diffusion  layer  (pumol  s-1) 

r 

reaction  rate  at  reaction  sites  (pmol  s-1) 

j 

diffusion  flux  (pimol  s-1  m-2) 

L a 

thickness  of  anode  diffusion  layer  (m) 

Lc 

thickness  of  cathode  diffusion  layer  (m) 

M 

mole  mass  (gmol-1) 

pb 

partial  pressure  in  gas  bulk  (atm) 

ptpb 

partial  pressure  in  the  immediate  vicinity  of 
tpbs  (atm) 

P 

pressure  (atm) 

pO 

standard  pressure  (atm) 

pO 

ru2o 

vapor  pressure  of  the  steam  at  the  temperature 
considered  (atm) 

R 

gas  constant  (=  82.05  x  10-5  Jmol-1  K-1) 

Ret 

charge  transfer  resistance  (Q) 

Rload 

load  resistance  (£2) 

Ro 

ohmic  resistance  (£2) 

S 

Laplace  operator 

T 

temperature  (K) 

V 

voltage  (V) 

fout 

fuel  cell  voltage  out  (V) 

Greek  letters 

£ 

porosity 

Odd 

activation  loss  (V) 

X 

tortuosity 

(E  vi) 

diffusional  volumes 

Superscripts 

0 

at  standard  pressure 

b 

gas  flow  bulk 

r 

flow  at  the  reaction  site 

s 

flow  at  the  upper  surface  of  diffusion  layer 

tpb 

triple  phase  boundary 

Subscripts 

12 

binary  diffusion 

a 

anode 

act 

activation 

c 

cathode 

ct 

charge  transfer 

eff 

effective 

h2 

hydrogen 

h2o 

water  vapor 

k 

Knudsen 

load 

load 

n2 

nitrogen 

0 

ohmic 

out 

output 

o2 

oxygen 

The  remainder  of  this  paper  is  organized  as  follows: 
the  principle  of  SOFC  is  discussed  in  Section  2.  The  dy¬ 
namic  modeling  is  discussed  in  Section  3.  Relational  pa¬ 
rameters  are  shown  in  Section  4.  Simulation  results  and 
analysis  are  given  in  Section  5,  followed  by  conclusions  in 
Section  6. 


2.  Fuel  cell  principles 

A  typical  H2-H20,Ni|YSZ|Pt,02  fuel  cell  is  shown  in 
Fig.  1.  The  chemical  reactions  inside  the  cell  that  are  directly 
involved  in  the  production  of  electricity  are 

at  anode  tpb:  H2  +  O2-  — >  H2O  +  2e-, 

at  cathode  tpb:  5O2  +  2e-  — >  O2-  (1) 


Hydrogen  Water  Heat  Electricity 


Fig.  1.  Principle  of  solid  oxide  fuel  cell. 
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2.1.  Voltage  output 

In  the  ideal  situation  (reversible),  electrical  work  is  equal 
to  the  Gibbs  free  energy  released  Agf,  i.e. 


Agf  =  —2FE 


The  electromotive  force  (EMF)  or  reversible  open  circuit 
voltage  E  of  the  hydrogen  fuel  cell  is  given  as 


Agf 
2  E 


where  F  is  Faraday  constant. 

Partial  pressures  or  gas  concentrations  affect  EMF  through 
the  Nernst  equation  [8]: 


E  =  E°  + 


RT  / «H2«of  \ 

\  flH20  J 


IF 


where  au2 ,  ao2  and  <2h2o  are  activities.  If  gases  behave  as 
ideal  gases,  the  activities  are 


P  h2  P  o2  Ph2o 

aH2  —  pQ  5  ^02  —  pQ  »  ^H20  —  Q 

r  r  ^h2o 


where  pn2,  po2,  and  pw2o  are  reactant  partial  pressures,  P° 
is  the  standard  pressure,  Pjj  Q  is  the  vapour  pressure  of  the 
steam  at  the  temperature  concerned.  E  is  also  called  Nernst 
voltage.  If  P°  =  Ph9o  =  h  the  equation  can  be  simplified  to 


E  =  E°  + 


RT 


In 


( 


1/2 

PU2  Po2 


2F  y  pu2o 


Irreversibilities  reduce  the  cells  voltage.  They  are  mainly  ac¬ 
tivation  loss,  ohmic  loss  and  concentration  loss.  The  voltage 
is  usually  modeled  in  the  steady  state  form  [8]: 


v  =  E-iRm  -A  In  Uj  -Blnfl-N  (6) 

where  i  is  the  current  produced  by  the  cell.  Rm  is  the  in¬ 
herent  resistance  of  the  fuel  cell,  z'o  is  the  exchange  current, 
an  important  parameter  of  weighting  the  activity  of  catalyst 
reaction.  i\  is  the  limiting  current,  at  which  the  fuel  is  used 
up  at  a  rate  equal  to  its  maximum  supply  rate.  A  and  B  are 
coefficients.  The  second  term  in  the  equation  represents  the 
ohmic  loss;  the  third  term  is  activation  loss  and  the  fourth 
term  represents  the  concentration  loss. 


The  maximum  current  that  a  cell  can  output  is  limited  by 
several  factors. 

The  first  one  is  reactant  supply  rates.  In  cell-level,  they  are 
controlled  by  concentration  gradients  between  tpbs  and  gas 
flow  bulks.  When  the  current  output  increases,  the  hydrogen 
and  oxygen  concentrations  at  tpbs  decrease  to  create  larger 
concentration  gradients.  Once  one  of  them  decreases  to  zero, 
the  supply  rate  reaches  its  maximum.  Current  can  not  be  in¬ 
creased  anymore,  and  it  is  the  maximum  current  fuel  cell  can 
provide.  Under  this  condition,  the  voltage  drops  to  zero  as 
described  by  Nernst  equation. 

Second,  current  output  is  also  limited  by  reaction  rates 
and  the  area  where  reactions  take  place.  In  most  cases,  reac¬ 
tions  are  fast  in  anode.  However,  because  cathode  reaction  is 
slow  [9],  current  output  is  also  limited  by  the  maximum  ion 
production  rate. 

Third,  current  output  is  controlled  by  voltage  and  load 
impedance. 

Fourth,  maximum  current  is  limited  by  ionic  conductivity 
of  electrolyte.  It  follows  Ohm’s  law.  The  limitation  resistance 
usually  merged  into  the  inherent  impedance.  In  normal  op¬ 
erating  ranges,  the  reactions  usually  do  not  reach  the  limits 
mentioned  before.  Thus,  the  current  is  determined  mainly  by 
Ohm’s  law. 

2.3.  Dynamic  process  of  fuel  cell 

When  fuel  gases  are  fed  to  cell,  they  must  flow  through 
the  gas  boundary  layers,  porous  support  layers  and  porous 
electrodes  to  tpbs,  where  reactions  take  place.  Neglecting 
the  reaction  dynamics,  these  mass  transport  processes  are  the 
main  dynamic  sources  of  the  fuel  cells  operated  under  the 
uniform  temperature  condition.  The  dynamics  of  fuel  cell 
including  rate  determination  steps  etc.  may  be  expressed  by 
transfer  function  block  diagram,  shown  in  Fig.  2. 

Another  important  source  that  affects  fuel  cell’s  dy¬ 
namic  behavior  is  its  inherent  impedance.  Fuel  cell  including 
SOFCs  inherent  impedance  is  well  studied  via  electrochemi¬ 
cal  impedance  spectroscopy  (EIS)  method  [10-16].  They  can 
be  simplified  as  equivalent  RC  circuits.  So  the  effect  of  in¬ 
herent  impedance  to  voltage  and  current  dynamic  behavior  is 
similar  to  a  RC  filter. 


3.  Modeling  approach 

3.1.  Assumptions 


2.2.  Current  output 

Neglecting  the  transit  dynamics  of  the  reactions,  the  rela¬ 
tion  between  current  and  reactions  can  be  expressed  as  [8]: 

*'  =  2FJ\h  =  2FJ\ho  =  4FJq2  (7) 

where  superscript  r  represents  the  fuel  consumption  or  water 
vapor  production  rate  at  tpbs. 


In  this  paper,  we  focus  on  dynamic  behavior  of  SOFC  due 
to  the  diffusion  and  impedance.  We  consider  two  scenarios: 
( 1)  an  elemental  volume  of  the  cell,  and  (2)  a  cell  with  uniform 
temperature.  The  results  derived  from  this  paper  are  valid  for 
both  scenarios.  This  serves  for  two  purposes:  (1)  when  con¬ 
sider  an  elemental  volume  of  SOFC,  the  assumption  of  uni¬ 
form  temperature  condition  is  valid  and  dynamic  behavior 
of  an  elemental  SOFC  can  be  used  as  a  building  block  for 
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Fig.  2.  Schematic  block  diagram  of  fuel  cell. 


a  complete  cell  dynamic  model  when  the  temperature  varies 
along  the  cell;  (2)  if  consider  a  whole  cell  operating  under 
uniform  temperature  instead,  we  actually  isolate  and  inves¬ 
tigate  SOFC  cell  dynamics  due  to  diffusion  and  impedance 
alone.  The  dynamic  responses  are  due  to  the  effect  of  diffu¬ 
sion  and  impedance  only  and  should  not  be  interpreted  as  the 
complete  dynamic  responses  of  the  cell. 

The  following  are  the  assumptions  if  we  consider  an  ele¬ 
mental  volume  of  the  cell: 

(1)  An  elemental  volume  of  H2-H20,Ni|YSZ|LSM,  air  cell 
[12].  The  definition  of  the  elemental  volume  in  different 
SOFC  cells  is  shown  in  Fig.  3. 

(2)  Cell’s  external  load  is  a  pure  resistance. 

(3)  Temperature  is  uniform  throughout  the  elemental  vol¬ 
ume. 

(4)  Gas  partial  pressures  in  flow  bulks  surrounding  the  ele¬ 
mental  volume  are  uniform. 


3.2.  I/O  variables 

The  aim  of  the  model  to  be  developed  is  to  describe  the  ex¬ 
ternal  characteristics.  So  input  variables  are  related  to  inputs 
of  SOFC.  Output  variables  are  cell’s  power  outputs,  such  as 
voltage,  current,  etc. 

The  input-output  variables  are  shown  in  Table  1 .  The  first 
input  variable  is  the  external  load.  Under  the  normal  operat¬ 
ing  conditions,  it  is  external  load  that  determines  the  current 
output  i  and  so  affects  the  reaction.  Different  load  impedance 
affects  the  output  properties  in  different  way.  In  order  to  in¬ 
vestigate  the  basic  dynamic  behavior  of  the  current  output, 
load  impedance  is  assumed  to  be  a  pure  resistance  7?ioad  in 
the  model. 

Other  input  variables  are  partial  pressures  of  reactants 
in  gas  bulks.  They  affect  the  reactant  supplies  and  reac¬ 
tion  rates  directly.  For  a  single  cell  under  normal  operating 
condition,  reactant  concentrations  or  partial  pressures  affect 
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Fig.  3.  Definition  of  finite  volume  in  different  SOFC  cells. 


the  diffusion  processes  and  the  Gibbs  free  energy,  and  thus 
the  voltage.  Therefore,  other  input  variables  are  p ^?,  p\^n 0, 

and  pf . 

Output  variables  of  the  model  include  potential  difference 
exerted  on  the  load  resistance,  Vout,  current  flow  through  the 
load  resistance,  i,  H2  and  O2  consumption  rates,  7^,  7q2, 
and  water  production  rate  /^q. 


3.3.  Diffusion 


The  path  of  mass  transport  from  the  flow  bulk  to  the  re¬ 
action  site  involves  two  stages.  First,  from  the  flow  bulk  to 
the  cell  surface  layer.  Second,  through  the  porous  electrode 
to  reaction  sites  [17].  In  the  first  stage,  mass  flux  diffuses 
through  the  boundary  layer  to  the  cell  surface.  In  the  sec¬ 
ond  stage,  mass  flux  diffuses  inside  the  porous  electrode.  In 
all  these  two  stages,  diffusion  is  the  main  means  of  mass 
transport. 

One  of  the  most  widely  used  diffusion  models  is  Fick’s 
law: 


dC 

dx 


Table  1 


Input  and  output  variables 


Inputs 

^load 

Load  resistance 

Ph2 

Partial  pressure  of  hydrogen  in  anode  gas  bulk 

Po2 

Partial  pressure  of  oxygen  in  cathode  gas  bulk 

pn2o 

Partial  pressure  of  water  vapor  in  anode  gas  bulk 

Outputs 

Vout 

Voltage  output  of  SOFC  cell 

i 

Current  output  of  SOFC  cell 

P 

JH2 

Hydrogen  consumption  rate  of  SOFC  cell 

p 

J02 

Oxygen  consumption  rate  of  SOFC  cell 

p 

Jn2o 

Hydrogen  consumption  rate  of  SOFC  cell 

The  mass  transport  equation  can  be  written  as 


dC  _  d2C 
dt  dx 2 


where  C  is  the  mass  concentration,  D  the  diffusion  coefficient, 
A  the  diffusion  cross-sectional  area,  and  v  is  the  diffusion 
depth  as  defined  in  Fig.  1. 

Fick’s  law  shows  that  concentration  is  dependent  on  dif¬ 
fusion  thickness.  In  order  to  get  the  concentrations  at  tpbs, 
one  usually  applies  the  finite  element  method  [18-22].  How¬ 
ever,  by  means  of  Laplace  transform,  the  analytical  dynamic 
relations  can  be  expressed  in  the  form  of  transfer  functions 
without  involving  spacial  variables. 


3.3.1.  Developing  mass  transport  transfer  function 

Perform  Laplace’s  transform  in  Eq.  (9)  to  convert  the  par¬ 
tial  differential  equation  to  the  ordinary  differential  equation 
[23]: 


d2CO)  .s' 

-  Dcw  = 0 

Boundary  conditions  are 

r  dC(Y) 

f(s)  =  -D—^\x= 0, 

dx 

Solving  Eq.  (10)  yields: 


Cb(s)  =  C(s)\x=L 


C(s)(x)  = 


exp  (vTL)  +  exP 

CbW_/gexp(yxL) 


exp 


+ 


“P  (vTL)  +exp(-vTL) 


exp 


(10) 


(11) 


5x 


X 


D 

(12) 


Assume  that  the  gases  are  ideal  gases  and  the  flow  area  is  A. 
At  tpbs  where  v  =  0,  the  partial  pressure  in  the  vicinity  of 
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tpb  is 


ypb(s> 


-exp(-^FO  i  rt jr(s) 

exP  (1/sL)  +  exP  {-\/t)L)JDs  a 


+ 


(13) 


On  the  surface  of  the  diffusion  layer,  mass  flow  rate  is  P  = 
AD(dC/dx)\x=L .  The  dynamic  relation  is 


2 


exP(VsL)  +exp(- 


+ 


exP(Vs~L)  ~exP 

exP  (VfL)  +exp  (-\/sL) 


(14) 


3.3.2.  Approximating  mass  transport  transfer  function 
Taylor’s  expansions  of  exp  {^/{sjD)L) 

exp(— y/{sjD)L)  are 


exp 


( 


=  1  + 


IS  9 

- L2  + 

2D 


1  S3/2  3 

6  D3/2  L 


+ 


1  .S'2 

24  D2 


L4  +  O 


exp 


1  S  2  1  S3/2  3 

_ T  1 _ T  5 

2D  6  D3/2 


and 


+ 


1  .S'2 

24  D2 


P  +  o 


Substituting  them  into  Eqs.  (13)  and  (14),  and  neglecting 
higher  order  term  yield: 


/7tpb(s)  =  G]pf(s)  +  Gpppb(s), 

Js(s)  =  GuJr(s)  +  Gp}ph(s)  (15) 

where 


partial  pressures  in  the  vicinity  of  tpbs  and  gas  flows  at 
electrode  surface  can  be  determined  uniquely  by  the  be¬ 
haviors  of  gas  consumption  rate  and  bulk  pressure  with¬ 
out  relying  on  concentration  distribution  along  the  diffusion 
path. 

In  this  dynamics  description,  only  two  parameters  are  in¬ 
volved.  The  first  one  is  the  thickness  of  the  diffusion  layer 
L.  It  depends  on  flow  velocity,  and  can  be  calculated  accord¬ 
ing  to  fluid  mechanics.  The  second  parameter  is  diffusion 
coefficient  D.  It  can  be  calculated  from  correlation  equations 
[24]. 


3.3.3.  Diffusion  coefficient 

In  porous  materials,  the  effective  diffusion  coefficient  is 
adjusted  [24]: 

Dcff=-D  (16) 

r 

where  s  is  the  porosity,  r  the  tortuosity  of  porous  materials, 
and  Dis  the  total  diffusion  coefficient. 

Considering  the  Knudsen  diffusion,  the  total  diffusion  co¬ 
efficient  is  [24] : 


1  _  1  ,  1 

d~d5+p 


(17) 


where  D 12  is  the  binary  diffusion  coefficient,  and  D k  is  the 
Knudsen  diffusion  coefficient.  If  pores  are  large  enough, 
Knudsen  diffusion  can  be  neglected  [24] . 

The  binary  diffusion  coefficient  D 12  is  modeled  by  the 
Fuller’s  correlation  [24]: 


f[(E-0!,',  +  (E‘v);/3]2 


(18) 


where  T  is  the  temperature,  M\  and  M2  are  the  mole  mass 
of  gases  1  and  2,  (fff  vf)\  and  vf) 2  are  the  diffusional 
volumes  of  gases  1  and  2,  respectively,  and  P  is  the  total 
pressure.  Good  agreement  between  the  Fuller’s  correlation 
and  measurements  has  been  reported  in  [25]. 


GJp  = 


GPP  = 


GpJ  = 


-5  ~  Ml  5 1 

l  +  Lfs+  2  A  ’ 

1  ^  2 D*  ^  24 D2 


1 


1  4-  ^S  +  -^-S2' 
1  ^  2 D*  ^  24 D2 


Gjj  = 


1 


1  +  2 
1  ^  2Z)^24D2 


Ls 


A 


1  +  rt 

1  ^  2 D*  ^  24 D2 


/?tpb  is  the  partial  pressure  in  the  vicinity  of  tpbs,  p°  the 

partial  pressure  in  gas  bulks,  P  the  gas  flow  into  the  outer 
surface  of  porous  material,  P  the  gas  consumption  or  wa¬ 
ter  production  rate  at  tpbs,  L  the  layer  thickness,  A  the  cell 
area,  D  the  effective  diffusion  coefficient,  R  the  gas  con¬ 
stant,  and  T  is  the  temperature.  So  dynamic  behavior  of 


3.4.  Voltage 


As  described  before,  fuel  cell’s  voltage  output  is  affected 
by  gas  partial  pressures  and  is  reduced  by  concentration  loss, 
activation  loss  and  ohmic  loss.  The  dynamic  behavior  of  the 
voltage  is  also  affected  by  these  factors. 

In  fact,  because  reactions  take  place  at  tpbs,  it  is  partial 
pressures  in  the  vicinity  of  tpbs  that  affect  the  EMF.  The  more 
appropriate  expression  of  Nernst  equation  should  be 


(19) 


where  p^,  p ,  and  p^0  are  the  partial  pressures  in  the 
vicinity  of  tpbs. 
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The  H2/O2  consumptions  lead  to  their  concentration  re¬ 
ductions  in  the  vicinity  of  tpbs.  The  voltage  drop  caused  by 
this  kind  of  concentration  reductions  is  named  as  concentra¬ 
tion  loss  in  the  literature  and  is  usually  corrected  by  a  static 
concentration  loss  term. 

The  activation  loss  is  the  potential  consumed  to  overcome 
the  activation  energy  barrier.  It  is  normally  described  by  the 
Butler- Volmer  correlation  [20]: 


1  =  i()  <  exp 


exp 


-(1  -  P) 


act 

RT 


(20) 


where  ft  is  the  transfer  coefficient  and  z‘o  the  exchange  current. 
The  transfer  coefficient  is  usually  taken  as  0.5  for  the  fuel  cell 
application  [20]. 

When  ft  =  0.5,  the  activation  loss  can  be  solved  from  Eq. 

(20): 


7? act, a 


7? act, c 


2  RT  1 

- sinh 

nF 

2  RT 

- sinh 

nF 


(21) 


The  exchange  currents  depend  on  the  activation  energy,  tem¬ 
perature,  and  partial  pressures.  They  can  be  calculated  as 
[17]: 


where  £act?a,  Fact,c  are  anode  and  cathode  activation  energy, 
respectively. 

Compensating  the  activation  loss,  the  irreversible  voltage 
is 


E  =  E°  + 


RT 

2F 


In 


^7act,a  ^7act,c 


with  E°  =  1.273  -  2.7645  x  10"4T  [17], 


(23) 


Rn  R, 


eta  Rf>  R(ic  Rc 


T°I°T 


z  inherent 


HI 

Cdla 


Cd,c 


Fig.  4.  Equivalent  circuit  of  inherent  impedance. 

electrolyte,  Ra  and  Rc  represent  the  ohmic  resistances  of  an¬ 
ode  and  cathode,  Rcta  and  Rctc  represent  the  charge  trans¬ 
fer  resistances,  respectively,  and  Cdia  and  Cdic  represent  the 
charge  double  layer  capacitors  between  anode,  cathode  and 
electrolyte. 

3.6.  Equivalent  circuit 

Theoretically,  impedance  spectra  plane  plot  indicates  two 
semicircles  for  the  impedance  shown  in  Fig.  4  [10].  In  SOFC, 
neglecting  diffusion  impedance,  these  two  semicircles  are 
merged  to  one  distorted  semicircle  [12].  This  means  that  the 
whole  inherent  impedance  can  be  approximated  by  one  RC 
unit.  In  another  aspect,  potential  difference  is  produced  be¬ 
tween  tpbs,  inside  the  two  metal  electrode  layers,  as  shown  in 
Fig.  1 .  That  means  voltage  fluctuation  due  to  potential  change 
is  smoothed  by  the  double  layer  capacitance.  So  a  reasonable 
equivalent  circuit  of  a  SOFC  cell  is  shown  in  Fig.  5,  where  R0 
is  the  total  ohmic  resistance  in  the  inherent  impedance,  Rct 
is  the  total  charge  transfer  resistance,  and  Cc t  is  the  approx¬ 
imated  charge  transfer  capacitance.  These  three  parameters 
can  be  identified  from  impedance  spectral  plan  plot  [10]  or 
Bode  plot  [12]. 


3.5.  Inherent  impedance 

When  current  flows  through  the  inherent  impedance,  volt¬ 
age  is  reduced.  The  voltage  drop  is  called  Ohmic  loss.  The 
dynamic  characteristics  of  the  inherent  impedance  also  affect 
the  dynamic  voltage  behavior  of  the  fuel  cell. 

Inherent  impedance  of  SOFC  is  complex.  Basically,  it 
consists  of  two  charge  capacitance  processes  and  three  re¬ 
sistance  processes:  anode  charge  double  layer  capacitance, 
cathode  charge  double  layer  capacitance,  ohmic  resistance, 
grain  boundary  resistances  and  electrode  reaction  resistances. 
Because  their  dynamic  behaviors  are  similar  to  RC  cir¬ 
cuits,  the  inherent  impedance  is  usually  modeled  as  equiv¬ 
alent  RC  circuits.  A  typical  equivalent  circuit  to  inherent 
impedance  is  shown  in  Fig.  4.  Here  Re  is  the  resistance  of 


Fig.  6.  Comparison  of  simulated  and  experiment  V-I  plot. 
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This  equivalent  circuit  is  different  from  the  one  shown 
in  [8],  where  the  ideal  battery  is  connected  in  serial  with  a 
RC  pair.  That  equivalent  circuit  shows  derivative  effect  on 
battery  potential  change.  The  derivative  phenomenon  on  fuel 
cell  has  not  been  observed  and  reported  in  literature.  The 
model  shown  in  Fig.  5  is  an  integral  circuit.  It  smooths  the 
voltage  output,  and  appears  to  be  more  reasonable. 

The  dynamic  behavior  of  the  voltage  output  based  on  the 
equivalent  circuit  is  therefore  determined  by 


1 

- E- 

tfctQt 


Tout  —  Vet  tR  o> 


1 


RctCct 


Vet  — 


1 


l  = 


Cct 

Vet 


R0  T  ^load 


(24) 


Hydrogen  consumption  rate: 


_^2^h2  +^l^H2  +^3“jTF^H2 


RT 


and  partial  pressure  in  the  vicinity  of  anode  tpb: 

..tpb  ,  tpb  ,  .tpb  ,  RT 

Ph2  =  ~h  I  Pu2  -  hl-Pu2  -hA  —  JH2 


4  RT  ..  h 

_  17a~AJH2  +hlPHi 

where  hi  =  24  D\jL\,  h2  =  12  DH2/L^, 

24D^2/Llh4  =  24Dii2/Ll. 

•  Oxygen  consumption  rate: 

Yh  =  ~°1  J02  -  °2Jo2  +  0\Jb2  +  T  /'o2 


(25) 


(26) 


(27) 


3.7.  Model  in  the  state -space  form 

Converting  the  transfer  function  form  of  dynamic  relations 
shown  in  Eq.  (15)  to  differential  equation  form,  the  dynamic 
behaviors  can  be  shown  as: 


and  partial  pressure  in  the  vicinity  of  cathode  tpb: 

RT  r 

°^Jo2 


..tpb 

Po2  = 


tpb 

°Wo2 


.  tpb 

°2Po2 


4  RT 


Lr  A 


Jro2+oiph0: 


(28) 


Fig.  7.  Step  responses  of  SOFC,  when  load  resistance  changes. 


40 


Y.  Qi  et  al.  /  Journal  of  Power  Sources  150  (2005)  32-47 


where  o\  =  24  D20JL/,  o2  =  \2D0l/Ll ,  o3  = 

24 Df/Lj  04  =  24Do.IL/ 

•  Water  vapor  production  rate: 

A 

^h2o  =  _u;i^h2o  -  w2Jb2o  +  w^Jn2o  +  w3j^Ph2o 

(29) 

and  partial  pressure  in  the  vicinity  of  anode  tpb: 

..tpb  tpb  .tpb  RT 

Pn2o  =  -  w2Pu2o  -  ^4  — 420 

4  RT  K 

-  ~  ^~^h2o  +  wiPn2o  (20) 

where  w\  =  24  D\i0/L\,  w2  =  12DH2o/^a’ 

m  =  24D^20/L^,  W4  =  24Dn2o/Ll- 


First  order  derivative  of  an  input  variable  can  be  approxi¬ 
mated  by  [26] : 

^  (  1  -  1-!— -  )  U(s)  (31) 

V  F  +  1/ 

in  the  differential  equation  form: 

it  =  Ku  —  v,  v  =  K2u  —  Kv  (32) 

where  v  is  a  intermediate  variable.  K  is  the  approximation 
factor,  usually  greater  than  10. 

Define  the  input  vector  as 

u  —  [  ^load  Pii2  Po2  Pu2o  ]T  (22) 

output  vector  as 

y  =  [  Vout  i  ^h2  ^o2  ^h2o  ]T  (24) 

Introduce  intermediate  variables  vh2,  vo2,  vh7o>  vr  and  de¬ 
fine  states  as 


x  —  [  ^ct  ^h2  ^h2  vh2  ^o2  ^o2  VC>2 


p  p 

Ju2o  jh2o 


tpb  .  tpb  tpb  .  tpb  tpb 

va2o  Pn2  Pu2  Po2  Po2  Pn2o 


.  tpb 

Pu2o 


(35) 


Fig.  8.  Step  responses  of  SOFC,  when  hydrogen  pressure  changes. 
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Combine  Eqs.  (23)-(32),  the  state- space  model  is  given  by 


or  in  a  compact  form: 

x  =  /(x,  u),  y  =  g(x,  u)  (36) 

4.  Parameters 

This  state-space  model  involves  several  parameters.  Nu¬ 
merical  values  and  their  sources  are  listed  in  Table  2. 
These  are  default  parameters  for  our  model.  However, 
since  these  parameters  are  collected  from  a  variety  of  ref¬ 
erences,  in  the  validation  of  the  model  with  a  specific 
data  set,  some  of  the  default  parameters  have  to  be  ad¬ 
justed  to  the  published  specific  parameters  for  a  meaningful 
comparison. 

5.  Simulation  results 


Simulation  is  done  according  to  the  developed  state-space 
model.  Steady  state  and  transient  behaviors  of  SOFC  at  dif¬ 
ferent  inputs  and  disturbances  are  investigated  through  the 
simulation. 

5.7.  Steady  state  output  and  model  validation 

Dynamic  data  that  are  useful  for  the  validation  of  the  devel¬ 
oped  model  are  very  limited  in  the  literature.  We  will  validate 


State  equations: 

1  1  1  JCl 

x\  =  - E - x\ 


7?ctCct  Rct^ct  Crt  u\  T  R 


o 


*2  =  *3 

1  x\  A 

±3  =  —  h\X2  —  I12X3  +  h\ - 1-  I13 - (Ku2~  X4) 


2 F  u\  T  Rq  RT 


±4  =  KZU2  —  Kx  4 


x5  =  x6 


1  x\  A 

*6  =  -o  1*5  -  02*6  +  01  — - — —  +  03  —  (Ku3  -  X~i) 

4 F  u\  +  Rq  RT 


±7  —  K  U3  —  Kx 7 


vg  =  V9 


±9  =  —w  1X8  —  W2Xg  +  W\ 


A 


1 


■xi 


2 F  \u\-\-R 


o 


+  W3—(KU4  -  Xio) 

K1 


±10  =  KZU4  —  Kx  10 


Xu  =  X\2 


X\2  = 


—h\X\\-h2X\2-h4 

RT  4  1 


RT  1  vi 


A  La  IF 

+  h\U2 


A  2  F  u  1  T  Rq 

±1  X\ 


u\  +  Ro  (u\  +  Rq)2 


(Ku  1  —  x\j) 


X\3  =  ^14 


*14  = 


■01X13—02X14  —  04 

RT  4  1 


RT  1  vi 


A  FC4F 

+  01U3 


A  4  F  u  1  T  Rq 

±1  Vl 


u\  4~  R0  (u  1  +  7?0)2 


(Ku  1  -  vi7) 


*15  =  *16 


RT  1 


V16  =  —w  1V15  —  W  2V16  —  w  4 


Vl 


RT  4  1 


A  2 F  \U\  +  Rq 
vi 


A  La  2F 
T  \JO\U4 

2 

±17  =  K  u  1  —  Kxyj 

Output  equations: 


-*  +_ 

u\  (ui  +  Rq)2 


(Ku  1  —  vn) 


R 


y  1  =  *1 


o 


u\  +  7?0 


Vl 


±2  = 


Vl 


w  1  T  R 


o 


±3  =  *2 


±4  =  *5 


y  5  =  *8 


42 


Y.  Qi  et  al.  /  Journal  of  Power  Sources  150  (2005)  32-47 


Table  2 

Parameters  used  in  simulation 


Symbol 

Description 

Source 

La  =  1  mm 

Thickness  of  anode  diffusion 
layer 

V-I  plot  [27] 

Lc  =  1  mm 

Thickness  of  cathode 
diffusion  layer 

V-I  plot  [27] 

A  =  1  cm2 

Fuel  cell  effective  area 

[27] 

Ret  0.9  Q 

Charge  transfer  resistance 

EIS  test  [12] 

R0  0.1  Q 

Ohmic  resistance 

EIS  test  [12] 

C  =  300  [xF 

Charge  transfer  capacitance 

EIS  test  [12] 

£act,a  =  llOkJmor1 

Anode  activation  energy 

[17] 

£act,a  =  120kJmol_1 

Cathode  activation  energy 

[17] 

0 

II 

CO 

Porosity 

[24] 

T  =  4 

Tortuosity 

[24] 

(E  V'W  =  7-07 

Diffusional  volume 

[24] 

(E  vdH2o  =  F2.1 

Diffusional  volume 

[24] 

(E  vdo2  =  16.6 

Diffusional  volume 

[24] 

(E  v-')n2  =  17.9 

Diffusional  volume 

[24] 

T  1223 K 

Work  temperature 

[12] 

p =0.97  atm 

Input  partial  pressure 

[12] 

Pair  =  1  atm 

Cathode  pressure 

[12] 

p^0  =0.03  atm 

Input  partial  pressure 

[12] 

time  [sec] 


the  model  according  to  the  experiment  results  of  Wagner  et 
al.  [12]  and  Tsai  and  Barnett  [27]  respectively. 

Wagner  et  al.  [12]  investigated  the  effect  of  diffusion  of 
SOFC,  and  the  dynamic  properties  are  shown  by  EIS  Bode 
plot.  Given  3%  humidified  H2,  pure  O2  and  infinite  large  load 
resistance  R\0SLd,  at  1223  K,  the  simulated  steady  state  OCV 
is  1 1 17  mV.  This  is  in  a  good  agreement  with  measured  OCV 
of  1 1 14  mV  by  Wagner  et  al.  [12]  under  the  same  condition. 
Because  Wagner  et  al.  [12]  did  not  present  the  thickness  of 
diffusion  layers,  the  experiment  V-I  plot  can  not  be  used  to 
verify  the  dynamic  model. 

The  comparison  of  the  simulated  steady  state  V-I  outputs 
and  the  experiment  result  of  Tsai  and  Barnett  [27]  is  shown 
in  Fig.  6. 

Tsai  and  Barnett  [27]  did  not  give  the  thickness  of 
the  diffusion  layer  directly.  When  i  reaches  the  limit 
i\,  at  which  the  gas  supply  rate  is  maximum,  thickness 
of  the  diffusion  layer  can  be  calculated  by  the  relation 
[27]: 

h  =  4 Ff0l  =  4FAD02  (37) 


Fig.  9.  Step  responses  of  SOFC,  when  air  pressure  changes. 
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Using  Eq.  (37),  the  thickness  is  calculated  as  Lc  =  1 .26  x 
10~3  m.  Because  anode  diffusion  does  not  contribute  to  the 
limitation,  the  thickness  in  anode  can  not  be  calculated.  La  is 
assumed  same  as  Lc  in  the  simulation.  The  inherent  resistance 
read  from  the  EIS  test  of  Tsai  and  Barnett  [27]  is  R0  =  0. 16  £2 
and  Rci  =  0.16  Q  at  1123  K. 

The  V-I  plot  shows  that  the  simulated  result  is  overall  in 
agreement  with  the  experiment  data  with  some  error.  A  main 
reason  for  the  error  is  the  unavoidable  model  parameter  er¬ 
ror  used  in  the  simulation.  They  are  collected  from  different 
sources  and  are  dispersed  in  a  wide  range  in  different  litera¬ 
tures. 

5.2.  Step  responses  due  to  R\0&&  step  changes 

As  shown  in  Fig.  7a,  when  load  resistance  has  step 
changes,  current  will  change  immediately  according  to 
Ohm’s  law.  Voltage  on  the  pure  ohmic  inherent  resistance 
R0  also  changes  at  once.  The  result  is  an  immediate  rise  of 
Vout-  Then,  for  the  reason  of  charge  transfer  capacity,  volt¬ 
age  on  the  charge  transfer  resistance  Rci  changes  slowly.  V0ut 


then  rises  slowly  to  the  final  value  [8].  The  response  can 
therefore  be  divided  into  three  stages,  immediate  change, 
then  fast  response  and  slow  response  following  the  fast 
response. 

For  the  specific  SOFC  investigated  by  Wagner  et  al.  [12], 
the  charge  transfer  capacity  is  around  several  hundred  micro¬ 
farads.  The  charge  transfer  resistance  is  around  1  Q.  So  the 
time  constant  is  less  than  1  ms.  The  effect  of  charge  transfer 
capacity  can  not  be  distinguished  from  the  immediate  rise  of 
voltage.  The  fast  response  in  Fig.  7a  is  actually  the  summation 
of  the  changes  of  voltage  drop  on  R0  and  Rct. 

The  slow  response  following  the  fast  response  is  caused 
by  the  slow  changes  of  concentration  at  the  vicinity  of  tpbs. 
When  Jx  changes  simultaneously  with  current  change,  j^tpb 
changes  at  a  slow  rate.  The  behaviors  are  determined  by 
Eqs.  (26),  (28)  and  (30).  Thus,  according  to  Nernst  re¬ 
lation,  E  responses  slowly  to  the  new  value.  Assume  the 
diffusion  thickness  La  =  Lc  =  1  mm,  simulated  time  con¬ 
stant  of  the  slower  voltage  response  is  around  0.02  s.  It  is 
in  good  agreement  with  Wagner’s  experiment  results  [12]. 
In  [12],  a  concentration  impedance  from  10  Hz  is  observed 
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Fig.  10.  Step  responses  of  SOFC,  when  water  vapor  pressure  in  anode  gas  changes. 
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in  the  EIS  test  Bode  plot.  The  time  constant  is  identified 
as  0.018  s. 

The  voltage  output  step  response  is  similar  to  the  cur¬ 
rent  interrupt  test  result  shown  in  [8],  where  the  dynamic 
behavior  of  Vout  is  explained  only  by  the  charge  transfer 
resistance  and  capacitance.  Contributions  from  concentra¬ 
tion  changes  in  the  vicinity  of  tpbs  are  not  noticed.  The 
simulation  results  show  that  it  is  concentration  changes  in 
the  vicinity  of  tpbs  that  are  the  main  reason  for  the  slow 
voltage  rise  behavior.  The  current  response  is  shown  in 
Fig.  7b  which  is  related  to  the  voltage  response  by  Ohm’s 
law. 

Reactant  consumption  rates  at  diffusion  layer  surface  re¬ 
spond  at  a  slow  rate,  as  shown  in  Fig.  7c-e.  Time  constant 
of  oxygen  is  different  from  that  of  hydrogen  and  water.  This 
different  response  processes  lead  to  different  behaviors  of  re¬ 
actant  partial  pressures  in  a  fuel  cell  stack  and  in  turn  affect 
the  electrical  power  output.  The  different  transient  responses 
will  also  lead  to  the  change  of  pressure  difference  between 
anode  side  and  cathode  side.  The  reason  for  these  slow  dy¬ 


namic  responses  is  the  existence  of  the  diffusion  layers.  The 
diffusion  layer  plays  a  role  like  a  buffer,  which  damps  the 
response  of  fuel  transportation  to  tpbs. 

5.3.  Reactant  partial  pressure  disturbances 

Partial  pressure  disturbances  in  gas  bulk  lead  to  electrical 
power  output  fluctuation  and  change  the  reaction  rate.  Step 
responses  due  to  hydrogen,  air,  and  water  vapor  pressures  are 
shown  in  Figs.  8-10,  respectively.  Simulation  results  show 
that  influences  from  pressure  disturbances  in  gas  bulks  on 
voltage  and  current  are  relatively  small.  That  is  because  the 
electrical  power  output  is  controlled  mainly  by  reaction  and 
load  resistance,  and  the  existence  of  the  buffer  effect  of  the 
diffusion  layer. 

5.4.  Effect  of  concentration  loss 

When  load  resistance  moves  to  a  small  value,  it 
leads  SOFC  to  the  concentration  loss  dominant  region, 


Fig.  11.  Step  responses  of  SOFC,  when  it  enters  concentration  loss  range. 
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the  dynamic  behaviors  of  SOFC  is  considerably  differ¬ 
ent  from  that  in  ohmic  loss  range,  as  shown  in  Fig. 
11.  The  transient  response  processes  behave  a  high  order 
dynamics. 

Simulation  result  shows  that  in  the  air  fueled  SOFC, 
it  is  the  oxygen  supply  that  is  the  main  resource  of 
the  concentration  loss.  Because  oxygen  fraction  in  air  is 
only  around  21%.  Oxygen  supply  also  plays  the  main 
role  in  the  dynamic  behavior  of  SOFC,  since  its  dif¬ 
fusion  coefficient  is  smaller  than  that  of  hydrogen  and 
water. 

5.5.  Effect  of  diffusion  layers 

Simulation  result  shows  that  the  thickness  of  the  diffu¬ 
sion  layer  has  strong  effect  on  properties  of  SOFC.  With  the 
increasing  of  the  thickness,  time  constants  of  the  dynamics 
increase.  The  resistance  from  diffusion  also  increases.  The 
effect  is  shown  in  Fig.  12. 


Diffusion  layers  in  fuel  cell  consist  of  not  only  electrode 
layer,  but  also  boundary  layer.  The  thicknesses  are  affected 
by  the  status  of  flow  bulks,  especially  the  flow  velocity.  So 
the  thicknesses  of  diffusion  layers  may  vary  within  a  large 
range. 

5.(5.  Effect  of  temperature 

Given  different  temperatures,  the  simulation  shows  that 
compared  to  its  effect  on  SOFCs  dynamics,  tempera¬ 
ture  has  larger  effect  on  the  voltage  output,  as  shown  in 
Fig.  13. 

Temperature  affects  the  Gibb’s  free  energy  Agf,  and  thus 
Nernst  voltage.  It  also  affects  the  conductivity  of  electrolyte, 
electrodes  and  connectors.  Simulation  shows  that  the  effect 
of  temperature  on  Gibb’s  free  energy  and  the  Nernst  voltage 
is  much  larger  than  on  other  factors. 

As  shown  in  Fig.  13,  the  effect  of  temperature  on  the  tran¬ 
sient  property  of  the  diffusion  response  is  negligible. 


Fig.  12.  Effect  of  diffusion  layer  thickness  on  step  responses  of  SOFC. 
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Fig.  13.  Effect  of  temperature 

6.  Conclusion 

A  dynamic  model  of  solid  oxide  fuel  cell  (SOFC)  with 
a  focus  on  the  diffusion  process  at  cell-level  is  derived  in 
the  paper.  The  species  dynamics  is  built  in  the  form  of  the 
state-space  model.  Dynamic  properties  of  SOFC  are  shown 
through  simulations.  It  demonstrates  that  diffusion  processes 
in  porous  layers  play  an  important  role  in  the  dynamic  be¬ 
havior  of  SOFC.  They  affect  concentrations  in  the  vicinity  of 
triple  phase  boundary  (tpb),  and  thus  the  electrical  properties. 

Simulation  shows  that  it  is  the  dynamic  behavior  of 
partial  pressures  in  the  vicinity  of  tpbs  contribute  more  to 
the  slow  rise  of  voltage  in  the  step  response  test  and  current 
interrupt  experiment,  not  the  charge  transfer  capacitance. 
Given  different  diffusion  thicknesses,  it  is  found  that  the 
voltage  output  and  dynamic  behavior  of  gas  consumption 
rates  are  affected  greatly  by  the  thicknesses  of  the  diffusion 
layers.  Simulations  also  indicates  that  temperature  has  large 
effect  on  the  voltage  output. 
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